Vacancy-stabilized crystalline order in hard cubes 
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We examine the effect of vacancies on the phase behavior and structure of systems consisting 
of hard cubes using event-driven molecular dynamics and Monte Carlo simulations. We find a 
first-order phase transition between a fluid and a simple cubic crystal phase that is stabilized by a 
surprisingly large number of vacancies, reaching a net vacancy concentration of ~ 6.4% near bulk 
coexistence. Remarkably, we find that vacancies increase the positional order in the system. Finally, 
we show that the vacancies are delocalized and therefore hard to detect. 
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The free energy of crystal phases is generally mini- 
mized by a finite fraction of point defects like vacan- 
cies and interstitials. Hovifcver, the equilibrium num- 
ber of such defects in most, if not all, colloidal and 
atomic/molecular crystals with a single constituent is 
extremely low. Exemplarily, for the face-centered cubic 
crystal of hard spheres, one of the few colloidal systems 
where the vacancy and interstitial fractions have been 
calculated, the fraction of vacancies and interstitials is on 
the order of 10^'* and 10^^ respectively, near coexistence 
[T]. As such, neither the vacancies nor the interstitials 
strongly affect the phase behaviour and so most studies 
of crystals ignore the effect of these defects. Nevertheless, 
vacancies and interstitials have a significant effect on the 
dynamics in an otherwise perfect crystal, as the main 
mechanism for particle diffusion is hopping of particles 
from filled to empty sites or between interstitial sites. 

In this paper, we examine a system of hard cubes 
where, as we will demonstrate, one cannot ignore the 
presence of vacancies. Arguably, a cube is one of the 
simplest non-spherical shapes and the archetype of a 
space-filling polyhedron. Surprisingly, despite the sim- 
plicity of this system, we find that the stable ordered 
phase is strongly affected by the presence of vacancies, 
so much that vacancies actually increase the positional 
order and change the melting point. Remarkably, the 
fraction of vacancies in this system is more than two 
orders of magnitude higher than that of hard spheres 
or any other known typical, experimentally realizable, 
single-component atomic or colloidal system, reaching 
6.4% near coexistence. Additionally, while purely hard 
(not rounded) colloidal cubes are yet to be realized, col- 
loidal cubes with various interactions are now a reality 
[2H7] and it is likely that hard cubes will be realized in 
the future. 

Here, we use Monte Carlo (MC) and event driven 
molecular dynamics (EDMD) [8j simulations to examine 
in detail the effect of vacancies on the equilibrium phase 
behavior of hard cubes. The model we study consists of 



N perfectly sharp hard cubes with edge length a in a vol- 
ume V. Aside from hard-core interactions which prevent 
configurations of overlapping cubes, the particles do not 
interact. In both types of simulation (MC and EDMD), 
overlaps are detected using an algorithm based on the 
separating axis theorem (e.g. Ref. [9 ). More informa- 
tion on the EDMD implementation for cubes is given in 
the Supporting Information. 

RESULTS 
Spontaneous vacancy formation 

The equation of state for a vacancy-free system of hard 
cubes has been the subject of a number of studies pHl - 
IT2] and was most recently examined by Agarwal and Es- 
cobedo [To]. It clearly shows a single, first-order phase 
transition between a fluid and an ordered phase. The 
authors of Ref. [TU] identified the ordered phase at co- 
existence to be a liquid crystal mesophase, i.e. a cubatic 
phase, which is characterized by the presence of long- 
range orientational order along three perpendicular axes, 
but a lack of long-range positional order. However, the 
authors noted that finite-size effects made it difficult to 
determine the extent of the positional order in their sys- 
tem, and based this identification on their observation of 
finite diffusion in the ordered phase. At high densities, 
both Refs. [10] and [12] agree that the ordered phase is 
a simple cubic crystal. 

There is no fundamental reason why diffusion cannot 
occur in a crystal, hence this is an insufficient criterion 
for distinguishing between a cubatic phase and a crys- 
tal. To study more closely the range of the positional 
order along the ordered branch we re-examined the in- 
termediate density region (near coexistence) using highly 
efficient EDMD [8.^ simulations allowing us to access sys- 
tem sizes more than an order of magnitude larger than 
the ones considered in Ref. [TU] • We simulated systems of 
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N = 40^ = 64000 particles starting from a simple cubic 
crystal lattice. Coexistence between the fluid and or- 
dered phase is observed directly for overall packing frac- 
tions 0.455 < 77 = Na^/V < 0.480. Snapshots of typical 
configurations are shown in the supporting information. 

Looking in detail at the EDMD simulations for rj be- 
tween 0.52 and 0.56, i.e., in the region where Ref. [TU] 
reported the cubatic phase, we noticed that in many sim- 
ulations the crystal lattice spontaneously transformed in 
one of two distinct ways (see the supporting informa- 
tion). In most cases, the simple cubic crystal did not 
maintain its original orientation in the box; instead it ro- 
tated, introducing defects and frustrations to the crystal 
lattice. In others, the system spontaneously added extra 
layers, i.e. extra lattice sites. As long as the crystal lat- 
tice remains aligned with the simulation box, the number 
of lattice sites can easily be measured from the number 
of peaks in the three-dimensional density profile of the 
cubes. Figure [l] shows a two-dimensional projection of 
such a density profile from simulations with N = 40'^ 
particles. From this plot we can determine that the sys- 
tem has Nl = 41"^ lattice sites. Thus, the system sponta- 
neously incorporated a large number of excess lattice sites 
into the crystal. The resulting crystal has a net vacancy 
concentration, of approximately a = {Nl—N)/Nl = 8%. 
It should be noted that since the volume and the number 
of particles in the system are fixed, it is generally not pos- 
sible to reach an equilibrium concentration of defects in 
the system. However, the formation of extra layers and 
lattice distortions both significantly increase the number 
of lattice sites in the crystal, and suggest that the ther- 
modynamically stable phase in this regime might be a 
vacancy-rich crystal structure. We note here that the 
systems sizes examined by Ref. pUl] would not have al- 
lowed for extra layers to form. Hence, we believe that the 
rotated, defective crystals we see are what the authors of 
Ref. [TD] identified as cubatic. 

To examine the effect of vacancies on the crystal struc- 
ture, we performed additional EDMD simulations on sys- 
tems with various net vacancy concentrations a — {Nl — 
N)/Nl for a number of packing fractions -q = Na'^ /V . 
The average global positional order in the system was 
measured using the positional order parameter averaged 
over all particles: 



{Gglobal} — 



(l/7V)^exp(zK.r,) 



(1) 



where K is a reciprocal lattice vector of the crystal un- 
der consideration and is the position of particle j. In 
all our plots we have chosen K to correspond to a single 
lattice vector, i.e Kg = with e = x,y,z and a the 
lattice spacing. To set the net vacancy concentration we 
placed iV = (1 — a)N]^ particles randomly on a simple 
cubic lattice with = 40'^ = 64000 lattice sites, and 
then rescaled the volume (and thus the lattice spacing) 



FIG. 1: Peaks in the density profile associated with a two- 
dimensional projection of the centers of mass of the cubes at 
packing fraction rj = 0.52, as measured in an EDMD simula- 
tion initialized with A'^ =40'' particles on a Nl =40^ simple 
cubic lattice, i.e. no vacancies. The number of lattice sites 
in the system spontaneously increased to Nl = 41'^ lattice 
sites, corresponding to a net vacancy concentration of ap- 
proximately 8%. The lines were added to highlight the 41 
evenly spaced layers in both (x and y) directions. 




FIG. 2: Global positional order (Ggiobai) as a function of the 
net vacancy concentration a for packing fractions rj — 0.51, 
0.52, 0.53, 0.54, and 0.56 in a system with Nl = 64000 lattice 
sites. The points indicate measurements of (Ggiobai) along the 
X, y, and z axes separately, while the lines are averaged over 
all three directions. 



of the box to the desired packing fraction 77. Hence, the 
resulting system has a lattice spacing that depends on the 
chosen packing fraction rj and net vacancy concentration 
a. These simulations show that there is a maximum in 
the global positional order as a function of net vacancy 
concentration a for varying packing fractions r] (Fig. [2|. 
An increase in order due to an increase in number of va- 
cancies is unexpected as typically the presence of defects 
(such as vacancies) decreases the order. This observa- 
tion suggests that adding defects reduces frustration in 
the crystal, potentially stabilizing a defect-rich crystal. 



3 



Defect realization 

The 'net vacancy concentration' a, defined above, is 
simply the excess of lattice sites Nl compared to the 
number of particles iV, divided by -/V^, i.e. the fraction 
of lattice sites that does not have a particle associated 
with it. In a typical system, for instance hard spheres, a 
vacancy is localized to a single lattice site and one can de- 
termine the number of vacancies by counting the number 
of empty lattice sites. In a hard-sphere crystal, a particle 
next to an empty lattice site is kept in place by its other 
neighbors. This is not the case for hard cubes and, as a 
result, the way vacancies manifest in this system is very 
atypical. In hard cubes, entirely empty lattice sites are 
rarely seen even in the vacancy-rich crystals near coexis- 
tence. Instead, a defect manifests itself as a finite-length 
chain of particles along one of the three principal axes in 
the crystal, in which the particles have a slightly larger 
inter-particle spacing than the average as shown in Fig. 
|4j Hence, if a vacancy extends over 4 lattice sites, as is 
the case for one of the vacancies highlighted in Fig. [4j 
then the vacancy is realized by 3 particles sharing 4 lat- 
tice sites in a one-dimensional chain. Additionally, while 
a two dimensional layer in a typical snapshot, such as 
in Fig. [4] shows regions of disorder, it should be noted 
that even at high vacancy concentrations, the crystal still 
shows a well-defined lattice spacing on average, which 
can be easily determined from the position of the peaks 
in the scattering function 5'(k) or the three-dimensional 
pair correlation function ^(r) (Fig. [3]). 

It should be noted that the net vacancy concentra- 
tion includes both vacancies as well as interstitials in the 
sense that each interstitial cancels a vacancy. However, 
the number of vacancies is higher than the number of 
interstitials resulting in the large positive net vacancy 
concentrations we find in this system. Similar to a va- 
cancy, interstitials are also not localized in this system 
and occur by n particles sharing n — 1 lattice sites. 

Phase diagram of hard cubes 

So far we have established a relationship between the 
order in the system and the net vacancy concentration. 
However, there is no way to determine the equilibrium 
net vacancy concentration and the phase boundaries from 
the EDMD simulations. A proper determination of the 
equilibrium concentration of vacancies as well as the 
phase diagram requires free-energy calculations. The free 
energy per particle (/ = /SF/N) of the solid with vacan- 
cies is given by: 

ntnW = /cin ( A) + Aot (A) + /co„,b , (2) 

where the first term is the translational free energy of 
a normal Einstein crystal |13j , the second term is the 
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FIG. 3: Three-dimensional pair correlation function g{x, y, z) 
with z = 0, measured in an EDMD simulation of a system 
of A'' = 64000 particles with packing fraction ri = 0.51 and 
vacancy concentration a = 0.055. The x, y, and z directions 
are taken along the three axes of the simulation box. The data 
is averaged over 50 snapshots, and over the four quadrants of 
the xy-plane. 

rotational free energy of the crystal [Hj, and the third 
term is the combinatorial entropy associated with placing 
N particles on Nl lattice sites: 

/comb = -{l/N) l0g[NLl/{Nl{NL - N)l)]. (3) 

Detailed descriptions on how these terms were calculated 
using Monte Carlo simulations are given in the Methods 
section. 

We determined the free energy as a function of net 
vacancy concentration a, at a fixed packing fraction 
T] — 0.52. As shown in the inset in Fig. [5|d, the min- 
imum in the free energy occurs for a high concentration 
of vacancies. Specifically, we find that at this density the 
number of particles is 4% lower than the number of lat- 
tice sites. While calculating the free energy as a function 
of a, we also observed that the free energy, excluding the 
combinatorial contribution, was almost linear. This is 
shown in the inset of Fig. [SJd. While we are uncertain to 
the origin of this linearity, it seems to suggest that the 
vacancies are only weakly interacting. 

For each value of a, the free energy as a function of 
density was obtained by combining the reference free en- 
ergies shown in Fig. [5]d with a separate equation of state 
measured for that value of a. By minimizing the resulting 
free energy with respect to a, we find that the number of 
excess lattice sites decreases as a function of the density, 
as expected (Fig. [sj:). 

Using a common tangent construction (see the sup- 
porting information) in combination with the determined 
free energies, we mapped out the phase diagram which 
is shown in Fig. [5^. We find coexistence between a fluid 
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FIG. 4: A typical configuration where three delocalized de- 
fects, and the particles directly surrounding the defects, have 
been highlighted (yellow) in a system with Nl = 8000 lat- 
tice sites at packing fraction rj = 0.56 and defect concentra- 
tion a = 0.016 from an EDMD simulation. The color of the 
other particles indicates the orientation of each particle with 
respect to the crystal lattice, with colors closer to blue indi- 
cating more closely aligned particles. In the defect furthest 
to the right, the highlighted area shows three cubes sharing 
four lattice sites. The uppermost defect has six cubes shar- 
ing seven lattice sites, and the bottom most defect has seven 
cubes sharing eight lattice sites. 



phase with coexisting density r]f = Na^/V = 0.45 and 
a vacancy-rich simple cubic crystal structure with coex- 
isting density rjm = 0.50 and net vacancy concentration 
a = (iVi — N)/N]^ ~ 0.064. The pressure and chemical 
potential at coexistence are f3pa^ = 6.16 and = 18.4, 
respectively. The inset of Fig. [Hj: shows the equations of 
state and phase transitions both including and excluding 
the effects of vacancies in the crystal phase. The presence 
of vacancies in the crystal significantly lowers the melting 
density compared to the one reported by Ref. 110| where 
they found that a defect-free crystal melted at rym = 0.52, 
while the freezing packing fraction is approximately the 
same. We note here that we also find a melting number 
density of 77™ = 0.52 if we exclude vacancies. Hence, we 
find that vacancies increase the range of stability of the 
simple cubic crystal. 



Diffusion 

As was already shown in Ref. [10], the ordered phase 
has appreciable diffusion in the intermediate density 
regime. This can be understood in terms of the delocal- 
ized defects, which diffuse through the crystal and allow 
particles to diffuse in the opposite direction. To investi- 
gate the effect of vacancies on the diffusion coefficient in 
the solid, we measured the long-time self-diffusion con- 
stant of cubes in the crystal phase using EDMD simula- 
tions. Figure [6] shows the diffusion constant as a function 
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FIG. 5: (a) The equilibrium phase diagram of hard cubes as 
a function of packing fraction r\. For r\ < 0.45 we predict a 
stable fiuid phase, for rj > 0.50 we find a stable crystal phase 
with vacancies (see, e.g. Fig. 2), and in between these two 
packing fractions we find coexistence between the crystal and 
fiuid. (b) Free energy per particle as a function of net vacancy 
concentration a at packing fraction rj = 0.52. The points 
correspond to measurements while the solid line is a guide 
to the eye. The estimated error in the free energies, based 
on independent runs, is ~ 0.004fc_BT'. Inset: Free energy per 
particle without taking into account the combinatorial free 
energy: fdd ~ f — ./comb- The solid line is a linear fit. (c) 
The net vacancy concentration a as a function of packing 
fraction r; in the crystal phase. The error bars are based 
on the width of the free energy minimum. Inset: Equations 
of state for the fluid phase (black), the stable vacancy-rich 
crystal (green), and the crystal without vacancies, i.e. a — 
(red, dashed). Note that the phase transition (dotted lines) 
shifts to lower densities when vacancies are taken into account. 
The phase transition for the vacancy-free system essentially 
coincides with the one in Ref. [10) . 



of density, where the net vacancy concentration at each 
density was chosen to correspond to the equilibrium net 
vacancy concentration shown in Fig. [5]:. Near coexis- 
tence, the diffusion coefficient increases significantly, up 
to a maximum of Dr/a^ = 0.05, where r = jSma'^ is 
the unit of time in the EDMD simulations and m is the 
mass of a single cube. For comparison, the diffusion con- 
stant in the fiuid at coexistence is Dr/a^ — 0.15, three 
times higher than in the coexisting solid. 

At fixed density, the diffusion constant increases ap- 
proximately linearly with the number of vacancies, with 
very little diffusion remaining at a = 0. An example of 
this is shown in the inset of Fig. [6j for packing frac- 
tion T] = 0.56. Note that even for vanishing net vacancy 
concentration, diffusion is still possible via the sponta- 
neous formation of delocalized interstitial-vacancy pairs. 
However, at the equilibrium net vacancy concentration 
(a = 0.013), the diffusion coefficient is eight times as 
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FIG. 6: Dimensionless diffusion coefBcient Dr/a^ in the solid 
phase as a function of the packing fraction rj. Here, r = 
•\/ (3ma^ is the unit of time in the EDMD simulations and m 
is the mass of a single cube. For each density, the crystal has 
the equilibrium net vacancy concentration as determined from 
the free energy calculations. The inset shows the diffusion 
coefficient at a fixed packing fraction rj = 0.56 and varying 
net vacancy concentration. 



high as in the vacancy-free crystal, indicating that va- 
cancies play a major role in the dynamics of the particles 
in the solid. 



DISCUSSION AND CONCLUSIONS 



similar result we are aware of is for parallel hard cubes; 
a somewhat artificial system which shows a very peculiar 
second order freezing transition from a fluid to a simple 
cubic crystal. [TBHT5] 

The existence of a vacancy-stabilized simple cubic 
phase in hard cubes leads to the question of whether 
vacancy-stabilized crystal structures are present in other 
anisotropic, entropy-driven systems. We would expect 
vacancies to be relevant for other (likely hard) systems 
with crystal structures where vacancies can delocalize. 
Note that delocalization also requires the absence of 
strong interactions that constrain particles to their lat- 
tice sites. For example, we would not expect high va- 
cancy concentrations to occur in the simple cubic struc- 
ture studied by Rechtsman et aL.|19j which resulted from 
isotropic interactions. Recently, the phase behavior of a 
large number of polyhedral shapes has been studied us- 
ing Monte Carlo simulations. [inillOlllI] Since vacancies 
are easily overlooked in the case of spontaneously formed 
crystals, and unlikely to form in simulations starting from 
a fully filled lattice, it is possible that high equilibrium 
vacancy concentrations occur in many of these systems. 
Specifically, for crystal structures where some or all of the 
neighboring particles can freely move into an empty lat- 
tice site, the possibility of crystal vacancies should likely 
be taken into account. Examples include the crystal 
structures predicted for hexagonal and triangular prisms 
in Ref. [lOj . or the 2d (rounded) hard squares studied in 
Refs. 



In this paper, we have examined the effects of vacan- 
cies on the phase diagram of hard cubes using both event 
driven molecular dynamics simulations as well as Monte 
Carlo simulations. From the molecular dynamics simu- 
lations it is clear that vacancies play an important role 
in the equilibrium phase behavior of hard cubes. Free- 
energy calculations show conclusively a first-order phase 
transition between a fluid phase and a vacancy-rich sim- 
ple cubic crystal phase with up to 6% vacancies. Up to 
the system sizes we have studied (40"^ particles), wc find 
long-range positional order for systems with an equilib- 
rium concentration of vacancies (see Fig. [s]). Thus, we 
find that the stable phase is a simple cubic crystal for all 
densities, albeit one with significant diffusion due to the 
high defect concentration. 

The number of vacancies in this system is orders of 
magnitude larger than typically seen in colloidal systems. 
The stability of this vacancy-rich phase can most likely 
be attributed to the delocalization of defects in the crys- 
tal (Fig. ^ : clearly, one vacancy provides additional free 
volume for multiple nearby particles, decreasing the en- 
tropic cost of creating a defect. In most other colloidal 
crystals, such as hard-sphere face-centered cubic crystals, 
particles near a vacancy are still confined to their lattice 
site by their remaining neighbors. As a result, the local 
entropy gain from a defect is much lower. The only other 



METHODS 
Event driven molecular dynamics simulations 

Please refer to the supporting text for a full descrip- 
tion. 



Free energy of the liquid 

Thermodynamic integration allows one to calculate the 
free energy for all densities assuming that both the equa- 
tion of state and the free energy at a reference density 
are known. When the free energy of a reference density 
F{pq) is known, the free energy as a function of num- 
ber density F(p) can be determined using the equation 
of state. In particular, the free energy is given by 

where p is the density and /3 = l/ksT with fc^ Boltz- 
mann's constant and T the temperature. To measure 
the free energy of the fiuid at a reference density, we 
used Widom insertion test particle method [TB]. The free 
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energy of the fluid at density po is then given by 

= PpiPo) (5 

Solid free energies with and without vacancies 

To calculate the Helmholtz free energy as a function 
of the density for the solid phase we use thermody- 
namic integration [13| in MC simulations of systems with 
Nl = 20^ = 8000 or iVi = 25 x 20 x 18 = 9000 lattice 
sites. We checked during our simulations that the num- 
ber of lattice sites did not change spontaneously. For 
systems with the same density and net vacancy concen- 
tration, the differences in free energy between these two 
lattice sizes were within the error of our measurements. 
However, while equivalent free energy calculations for a 
smaller system (A'"^ = 1000 lattice sites) yielded quali- 
tatively similar results, finite-size effects were noticeable 
when compared to the larger systems. 

For the reference free energy of a crystal without va- 
cancies, we use a variation on the method introduced 
by Frenkel and Ladd |T3], where particles are tied to 
their respective lattice sites with springs, transforming 
the crystal into a non-interacting Einstein crystal for a 
sufficiently high spring constant A. In this case, we also 
add an aligning potential to handle the orientational de- 
grees of freedom of the particles [T^. Using the same 
coupling constant A that attaches the particles to their 
lattice sites, the aligning potential is given by 

N 
i—l 

where x(y) is a unit vector along the x{y) axis, and u^j-, 
with j = 1,2,3, are three mutually perpendicular face 
normals associated with particle i. Also, f3 — l/ksT is 
the inverse thermal energy, where fc^ is Boltzmann's con- 
stant and T the temperature. The parameter A controls 
the strength of the external potentials; hence for A = 
the system reduces to pure hard cubes, and for A — Xm 
with Am sufficiently large, the particles in the crystal are 
non- interacting. 

To calculate the free energy of a system with vacancies, 
instead of fixing the particles to a specific lattice site, we 
attach the particles to their nearest lattice site [TS] using 

C/ext(A) = A^ ^ |r, - r"(r,)n + C/rot(A) (7) 

where r°(ri) is the position of the lattice site nearest 
to r^. In this case, the dimensionless free energy per 
particle (/ — /3F/N) of the noninteracting system is 
fctnW = /ein(A) -h/rot(A) -h/comb where the first term is 
the translational free energy of a normal Einstein crystal 



|13) . the second term is the rotational free energy of the 
crystal [M], and the third term is the combinatorial en- 
tropy associated with placing N particles on Nl lattice 
sites: Umb = -{l/N)log[NL\/iNl{NL - Ny.)]. The full 
free energy of the crystal of hard cubes with vacancies is 
then given by: 

In contrast to the free energy calculations for systems 
without vacancies [13 , the center of mass of the sys- 
tem is not fixed in these simulations. To equilibrate the 
position of the center of mass, we introduce MC moves 
that collectively translate every particle in the system 
|15j . Additionally, moves that translate a single particle 
by exactly one lattice vector are introduced in order to 
improve sampling of different distributions of vacancies 
over the crystal. For a system with full lattice site occu- 
pancy {N = Nl) and thus no vacancies, we obtain good 
agreement between the two methods. 

Diffusion constants 

To measure the long-time self-diffusion constant in the 
crystalline phase shown in Fig. [6] we performed EDMD 
simulations in systems of Nl — 8000 lattice sites, for a 
range of densities. The vacancy concentration was chosen 
to correspond to the equilibrium vacancy concentration 
shown in Fig. [s];. The diffusion constant was calculated 
from the slope of the mean squared displacement as a 
function of time. An example of a plot showing the mean 
squared displacement is shown in the supporting text. 
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